	function [GDcolat,GDlon] = GG2GMD(ggcolat, gglon)
% Transformation from Geographic Coordinates to Geomagnetic Dipole Coordinates. 
%*********************************************************
% The pole position was calculated as the average of 
% the north pole position calculated for years 2000 
% and 2001 using IGRF2000.
%*********************************************************
%  input:	GGcolat	= geographic colatitude
%		GGlon	= geographic longitude (0-360)
% output:	GDcolat	= geomagnetic dipole colatitude
%		GDlon	= geomagnetic dipole longitude (0-360)
%********The output is a rough estimate of GM coords ******
% Geographic coordinates: GEO
%     X = intersection of Greenwich meridian and geographic equator
%     Z = Geographic North Pole
% Geomagnetic Dipole coordinates: MAG
%	Y = Intersection between geographic equator and the geographic
%		meridian 90deg East of the meridian containing the dipole
%		axis
%	Z = Dipole axis
% the third axis completes the Right hand triad.
% Definitions from Hapgood 1992.
%--------------------------------------------------------------------
% It is not recommended to use gdlon in geomagnetic studies because
% most variations in the field are time dependent. Time is fixed to 
% Geographic longitude.
%   transformation is rm = M.rg
    
%			|cos(a)cos(b)	cos(a)sin(b)	-sin(a)	|
%       M =	|  -sin(b)		    cos(b)		  0		|
%			|sin(a)cos(b)	sin(a)sin(b)	cos(a)	|
% where the geographic coordinates of the North Pole are 
%				((90-a) N, b E)  a=NPcolat b=NPlon

%			| sin(c)cos(d)	|  
%	 rg =	| sin(c)sin(d)	|
%			|     cos(c)	|
%   c= GGcolat d = GGlon
%			| sin(e)cos(f)	|	
%	 rm =	| sin(e)sin(f)	|	tan(GDlon) = sin(f)/cos(f)
%			|     cos(e)	|   cos(GDlon) = cos(e)

% Hence........
% tan(GDLon) = cos(b)sin(c)sin(d) - sin(b)sin(c)cos(d)
%              ---------------------------------------
%   cos(a)cos(b)sin(c)cos(d)+ cos(a)sin(b)sin(c)sin(d) - sin(a)cos(c)

% cos(GDcolat) = sin(a)cos(b)sin(c)cos(d) + sin(a)sin(b)sin(c)sin(d) 
%                 + cos(a)cos(c)
% BEGIN PROGRAM
% radiuns conversion
rad = pi/180;
% colatitude of north pole igrf2000 (90.0 - 81.70)
	NPcolat = 8.3*rad;
% longitude of north pole igrf2000
	NPlon = 277.36*rad;
    GGcolat = ggcolat*rad;
    GGlon= gglon*rad;
% initialize variables
	dividend = [];
	divisor = [];
	GDlon = [];
	GDcolat = [];
% determine dividend of GDlon
	dividend = cos(NPlon).*sin(GGcolat).*sin(GGlon);
	dividend = dividend - ( sin(NPlon).*sin(GGcolat).*cos(GGlon) );
% determine divisor of GDlon
	divisor = cos(NPcolat).*cos(NPlon).*sin(GGcolat).*cos(GGlon);
	divisor = divisor + ...
               ( cos(NPcolat).*sin(NPlon).*sin(GGcolat).*sin(GGlon) ) ;
	divisor = divisor - ( sin(NPcolat).*cos(GGcolat) );
% determine GDlon
	GDlon = atan(dividend./divisor);
% determine GDcolat
	GDcolat = sin(NPcolat).*cos(NPlon).*sin(GGcolat).*cos(GGlon);
	GDcolat = GDcolat + ...
            ( sin(NPcolat).*sin(NPlon).*sin(GGcolat).*sin(GGlon) );
	GDcolat = GDcolat + ( cos(NPcolat).*cos(GGcolat) );
	GDcolat = acos(GDcolat);
% convert from radiuns to degrees
    GDcolat = GDcolat/rad;
    GDlon = GDlon/rad;
%% END of FILE %%%